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Proton  exchange  membrane  fuel  cell  (PEMFC)  systems  with  their  own  fuel  conversion  unit  typically 
consist  of  a  fuel  processing  subsystem,  a  fuel  cell  stack  subsystem,  a  work  recovery-air  supply  subsys¬ 
tem,  and  a  power  electronics  subsystem.  Since  these  subsystems  have  different  physical  characteristics, 
their  integration  into  a  single  system  level  unit  makes  the  optimization  problem  of  synthesis/design  and 
operation/control  highly  complex.  Thus,  dynamic  system/subsystem/component  modeling  and  highly 
effective  optimization  strategies  are  required.  Furthermore,  uncertainties  in  the  results  of  system  opti¬ 
mization  can  be  affected  by  any  number  of  sources  of  uncertainty  such  as  the  load  profiles  and  cost 
models.  These  uncertainties  can  be  taken  into  account  by  treating  the  problem  probabilistically.  The 
difficulty  with  doing  this,  particularly  when  large-scale  dynamic  optimization  with  a  large  number  of 
degrees  of  freedom  is  being  used,  is  that  the  traditional  probabilistic  approaches  are  simply  too  com¬ 
putationally  intensive.  This  difficulty  can  be  overcome  by  the  use  of  approximate  approaches  such  as 
the  response  sensitivity  analysis  (RSA)  method  based  on  Taylor  series  expansion.  Thus,  in  this  paper,  a 
stochastic  modeling  and  uncertainty  analysis  methodology  for  energy  system  synthesis/design  and  oper¬ 
ation/control  which  uses  the  RSA  method  is  proposed  and  employed  for  calculating  the  uncertainties  on 
the  system  outputs.  Their  effects  on  the  development  and  control  optimization  of  a  5  kWe  PEMFC  system 
are  assessed  by  taking  the  uncertainties  into  account  in  the  objectives  and  constraints.  It  is  shown  that 
these  uncertainties  significantly  affect  the  reliability  of  being  able  to  meet  certain  constraints  (e.g.,  that 
on  the  CO  concentration)  during  the  synthesis/design  and  operation/control  optimization  process. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Fuel  cells  are  promising  candidates  as  alternative  energy  con¬ 
version  devices  for  transportation,  stationary  power,  and  portable 
applications.  Over  the  last  several  decades,  there  has  been  an  ever 
increasing  interest  worldwide  in  the  development  and  use  of  fuel 
cells  because  they  provide  an  efficient,  safe,  and  reliable  power 
solution.  However,  although  fuel  cell  systems  (FCSs)  exhibit  these 
great  benefits,  there  still  exist  a  number  of  technical  barriers  which 
must  be  overcome  such  as  material  durability,  reliable  interfacing 
with  conventional  utility  grids,  commercial  viability,  etc.  Further¬ 
more,  to  realize  such  systems  and  the  potential  benefits,  there  is  a 
need  for  optimal  system  synthesis/design  and  operation/control. 

The  optimal  synthesis/design  and  operation/control  of  FCSs 
requires  advanced  techniques  for  being  able  to  determine  the 
syntheses/designs  and  dynamic  operating  stages  of  such  systems, 
which  make  them  technically  and  economically  viable,  because 
FCSs  are  highly  non-linear  systems  that  have  a  number  of  reliability 
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issues  (e.g.,  catalyst  poisoning,  structural  degradation,  and  temper¬ 
ature  and  pressure  limitations),  which  must  be  addressed  in  order 
to  avoid  system  failures.  FCSs  typically  consist  of  a  fuel  processing 
subsystem  (FPS),  a  fuel  cell  stack  subsystem  (SS),  a  work  recovery- 
air  supply  subsystem  (WRAS),  and  a  power  electronics  subsystem 
(PES).  Only  the  first  three  of  these  are  considered  in  this  study. 
Since  each  subsystem  has  significantly  different  physical  charac¬ 
teristics  (i.e.  thermodynamic,  kinetic,  geometric,  temporal,  etc.), 
their  integration  into  a  single  optimal  system  renders  the  prob¬ 
lems  of  dynamic  system  synthesis/design  and  operation/control  to 
be  highly  complex.  Moreover,  it  is  also  very  important  to  synthe¬ 
size/design  and  control  the  systems  intelligently  in  order  not  only 
to  avoid,  as  mentioned  above,  episodes  of  system  failure  but  as 
well  to  obtain  optimal  system  operation  across  an  entire  load  spec¬ 
trum  both  in  terms  of  maximizing  energy  savings  and  minimizing 
environmental  effects. 

Conventional  approaches  for  synthesis/design  are  based  on  a 
single,  full  load  condition  at  steady-state.  This  can  significantly  sim¬ 
plify  the  system  optimization  problem  involved.  However,  because 
the  optimization  reflects  only  a  single  condition,  it  may  provide 
over-/underestimated  solutions  of  synthesis/design  and  also  lead 
to  non-optimal  solutions  for  operation/control.  Therefore,  with  the 
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Nomenclature 

C 

capital  cost  ($) 

Qz 

double  layer  capacitance  in  cathode  (Fcnrr2) 

Cp 

specific  heat  (kj  kmol-1  K-1 ) 

Enet 

system  net  power  generation  (W) 

Ewras 

motor  parasitic  power  by  WRAS  (W) 

t 

vector  of  inequality  constraints 

7? 

vector  of  equality  constraints 

m 

bulk  flow  proton  concentration  (mol  cm-2) 

[H+]  o 

reference  proton  concentration  (mol cm-2) 

j 

cell  current  density  (A cm-2) 

Jr 

reaction  current  density  (A cm-2) 

771 

mass  flow  rate  (kgs-1) 

n 

number  of  electrons 

77 

Molar  flow  rate  (mol  s-1 ) 

P 

pressure  (bar) 

P 

power  (W) 

T 

temperature  (I<) 

wM 

required  load  on  motor  (W) 

WE 

work  recovered  by  expander  (W) 

Wc 

required  work  on  compressor  (W) 

Greek  letters 

a 

transfer  coefficient 

K 

specific  heat  ratio 

X 

coupling  function 

m 

system/component  efficiency  of  i 

p 

voltage  overpotential 

V 

variance 

8 

standard  deviation 

Superscripts 

design 

design  case 

in 

inlet 

out 

outlet 

exception  of  systems  which  always  operate  at  a  single  load  point,  it 
is  necessary  to  take  into  account  part  as  well  as  full  load  conditions 
during  the  development  phase  of  a  system. 

Another  difficulty  is  that  as  energy  systems  become  larger  and 
more  complex,  the  greater  number  of  possible  system  configu¬ 
rations  and  technologies  that  could  possibly  meet  the  designer’s 
objectives  optimally  increases  greatly.  In  addition,  the  system  may 
need  to  be  developed  taking  into  account  both  the  transient  and 
environmental  effects  on  system  performance.  Thus,  the  difficulty 
of  developing  the  entire  system  via  the  formulation  of  a  single 
optimization  problem  in  which  the  optimal  synthesis/design  and 
operation/control  of  the  system  are  achieved  simultaneously  is 
great  and  rather  problematic  due  to  the  complexities  involved. 
These  complexities  are  further  heightened  with  the  introduction 
of  uncertainty  quantification  into  the  problem,  transforming  the 
problem  from  a  purely  deterministic  one  into  a  probabilistic  one. 
Uncertainties,  system  complexity  and  non-linearity,  and  large 
numbers  of  decision  variables  quickly  render  the  single  optimiza¬ 
tion  problem  unsolvable  by  conventional  single-level  optimization 
strategies.  They  can,  nonetheless,  be  handled  by  sophisticated 
multi-level  optimization  strategies  (i.e.,  decomposition  strategies). 
Decomposition  breaks  the  large-scale  optimization  problem 
down  into  a  set  of  approximately  equivalent  smaller  optimization 
problems  in  order  to  facilitate  the  optimization  procedure.  Decom¬ 
position  approaches  are  very  effective  for  the  optimization  of 
dynamic  systems  which  have  highly  nonlinear  characteristics  with 
a  large  number  of  degrees  of  freedom.  In  this  research,  physical 


Fig.  1.  Schematic  of  general  modeling  and  optimization  under  uncertainty. 


decomposition  [1]  is  applied  to  the  dynamic  synthesis/design  and 
operation/control  optimization  problem  of  a  PEMFC  system. 

As  to  the  uncertainties  in  the  system  synthesis/design  and  oper¬ 
ation/control  optimization  results,  these  can  be  due  to  any  number 
of  sources,  which  can  be  categorized  into  direct  ones  such  as  those 
which  result  from  computational  errors  and  indirect  ones  such 
as  those  which  result  from  a  load  profile  and  the  fuel  cost  [2]. 
This  study  has  focused  on  evaluating  the  uncertainties  in  system 
responses  due  to  indirect  uncertainty  sources  because  the  sys¬ 
tem  optimization  is  significantly  influenced  by  the  load  profile  and 
cost  information.  Quantifying  these  uncertainties  then  becomes  an 
important  task  in  the  development  of  the  system.  In  this  study, 
response  sensitivity  analysis  (RSA)  is  employed  and  developed  for 
use  with  dynamic  energy  system  optimization.  The  load  profile,  fuel 
cost,  and  cost  models  are  treated  as  probabilistic  input  values  and 
uncertainties  in  the  output  results  are  determined. 

Fig.  1  describes  how  all  of  the  issues  of  system  modeling,  of 
multi-level  optimization,  of  the  choice  of  optimizer,  of  uncertainty 
quantification,  etc.  just  outlined  are  integrated  into  a  complete 
process,  which  simultaneously  determines  the  dynamic  system 
synthesis/design  and  operation/control  optimization.  The  purpose 
of  this  research  is  to  thoroughly  study  these  issues  and  develop  and 
apply  appropriate  techniques  to  address  them  in  the  development 
of  a  5  kWe  PEMFC  system. 


2.  System  description  and  modeling 

2A.  PEMFC  system  description 

Fig.  2  depicts  the  initial  non-optimized  configuration  of  a  5  kWe 
PEMFC  system  developed  here,  taking  into  account  all  of  the  equip¬ 
ment  and  recovery  loops  necessary  for  maximizing  total  system 
electrical  efficiency.  The  main  objective  of  the  FPS  is  to  provide  the 
hydrogen  rich  gas  required  for  the  operation  of  the  PEMFC  stack.  In 
particular,  the  FPS  consists  of  three  main  steps;  a  desulfurization 
step  as  a  fuel  (e.g.,  natural  gas)  preparation  process,  a  reforming 
step  to  reform  the  hydrocarbon  fuel  to  a  hydrogen  rich  gas,  and  a 
CO  removal  step  to  reduce  the  CO  concentration  level  to  less  than 
10  ppm.  Organic  sulfur  contained  in  natural  gas  must  be  removed 
by  the  desulfurizer  to  prevent  the  detrimental  effects  on  the  cata¬ 
lysts  in  a  typical  PEMFC  posed  by  the  presence  of  sulfur. 

For  the  reformer  model,  steam  reforming  technology  is  selected, 
and  a  steam  generator  and  a  combustor  are  introduced  to  sup¬ 
port  the  steam-methane  reformer  (SMR).  Reformate  gas  from  the 
SMR  usually  contains  between  an  8  to  10%  concentration  of  CO. 
It  is  necessary  to  reduce  the  CO  concentration  below  10  ppm  for 
safe  operation  of  the  PEMFC  system.  To  satisfy  this  requirement,  a 
three-step  CO  removal  unit  is  designed  for  this  FPS  which  includes 
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FPS:  fuel  processing  subsystem 

WRAS:  work-recovery  and  air-supply  subsystem 

SS:  stack  subsystem 

PES:  power  electronics  subsystem 


C:  compressor 
E: expander 
MIX:  mixing  chamber 
HX:  heat  exchanger 


SMR:  steam  methane  reformer 
HTSR:  high-temperature  water-gas  shift  reactor 
LTSR:  low-temperature  water-gas  shift  reactor 
PrOx:  preferential  oxidation 
P-HX:  pre-heat  exchanger 


Fig.  2.  Proposed  non-optimized  5  kWe  PEMFC  system  configuration. 


high-/low-temperature  water-gas  shift  reactors  (HTSR/LTSR)  and 
a  preferential  oxidation  (PrOx)  reactor. 

As  to  the  SS,  a  standard  PEMFC  stack  based  on  Du  Pont’s 
Nation™  is  used  which  requires  humidification,  temperatures  not 
in  excess  of  80  °C,  and  levels  of  CO  not  in  excess  of  10  ppm  to 
avoid  poisoning  of  the  electrochemical  catalyst.  Thermal  manage¬ 
ment  in  the  SS  is  addressed  by  integrating  the  PEMFC  stack  with  a 
water  cooling  cycle  that  controls  the  operating  temperature  of  the 
stack. 

The  WRAS  plays  a  significant  role  in  the  energy  integration  of  the 
whole  system.  It  consists  of  a  compressor,  an  expander,  and  a  motor. 
The  compressor  provides  compressed  air  to  the  fuel  cell  stack  and 
other  balance  of  plant  components  and  is  driven  by  an  expander 
and  motor.  An  electric  motor  is  used  to  supply  additional  power  to 
the  compressor  in  case  the  power  extracted  from  the  expander  is 
not  enough  to  run  the  compressor. 

2.2.  Subsystem  modeling 


The  current  that  is  created  by  these  electrochemical  reactions 
depends  directly  on  the  potential  difference  between  the  carbon 
support  for  the  catalyst  and  the  polymer  material  surrounding 
the  carbon  supported  catalyst  as  well  as  on  the  reactant  concen¬ 
trations.  The  reaction  current  density  jr  can  be  expressed  by  the 
Butler-Volmer  equation  as 


jr  =j  (A 


Pi  [H+] 

p?  [H+]0 


-  exp  - 


(-^ ) 


(1) 


where  [H+]  is  the  bulk  flow  proton  concentration,  [H+]0  the  ref¬ 
erence  proton  concentration,  rj  the  voltage  overpotential  due  to 
activation  and  concentration  losses,  a  the  transfer  coefficient,  n 
the  number  of  electrons  of  the  elementary  electrochemical  reac¬ 
tion,  Ar  the  effective  catalyst  area  per  unit  geometric  area,  and  j0 
the  exchange  current  density.  The  cell  current  density,  j,  is  the  sum 
of  the  reaction  current  density,  jr,  and  the  charge  storage  which  rep¬ 
resents  stored  charge  in  the  electrical  double  layer  at  the  interface 
between  the  catalyst  layer  and  the  diffusion  layer,  i.e., 


2.2.1.  SS 

A  semi-empirical,  one-dimensional  approach  based  on  Ceraolo 
et  al.  [3]  is  used  to  model  both  the  steady  state  and  transient 
behavior  of  a  single  PEMFC  as  well  as  fuel  cell  stack.  Each  cell 
utilizes  an  H2-rich  gas  mixture  as  the  fuel  and  air  as  the  oxi¬ 
dant,  both  humidified.  The  stack  size  (i.e.,  the  number  of  cells 
and  cell  dimensions)  is  determined  via  an  optimization  process, 
which  results  in  a  system  and  set  of  components  that  meet  the 
power  demands  and  design  objectives.  The  temperature  of  the 
fuel  cell  stack  is  assumed  to  be  uniform  (a  good  assumption  as 
the  temperature  variation  across  the  stack  is  relatively  small)  and 
is  used  to  determine  the  flow  rate  of  the  stack  cooling  water 
used  to  maintain  this  stack  temperature.  The  water  vapor  con¬ 
tained  in  the  reactant  mixtures  in  the  pores  of  the  cathode-side  gas 
diffusion  and  electrode-catalysts  layers  is  assumed  to  be  in  equi¬ 
librium  with  the  surrounding  liquid  phase  so  that,  consequently, 
the  partial  water  pressure  is  uniform  throughout  these  layers.  Fur¬ 
thermore,  the  membrane  electrolyte  is  assumed  to  be  saturated 
completely  with  water  so  that  its  conductivity  is  only  a  function  of 
temperature. 


•  •  r  A 

J~Jr  +  Cdldt 


(2) 


where  Q  b  the  double  layer  capacitance,  is  assumed  constant. 

Ceraolo  et  al.  [3]  propose  the  following  semi-empirical  expres¬ 
sion  for  the  exchange  current  density: 


joA 


r(? 


bh) 


X  exp i-a-ij-aij5) 


(3) 


With  Eqs.  (2)  and  (3),  transient  behaviors  of  the  cell  current 
density  can  be  calculated  as  a  function  of  time. 

2.2.2.  WRAS 

A  screw  type  compressor  and  expander  are  selected  for  the 
WRAS.  A  screw  compressor  can  provide  a  high  discharge  pressure 
at  low  mass  flow  rates  and  low  speeds  and  also  provides  the 
oil-free  output  required  by  fuel  cell  systems.  Therefore,  even 
though  screw  compressors  are  more  expensive  than  other  typical 
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Rotor  speed  factor  (rpm/T°  5) 


400  600  800  1000  1200 


Table  1 

Economic  assumptions  for  the  economic  analysis. 


Parameter 

Values 

N units 

Production  volume 

500,000 

h 

Operating  hours  per  year 

7920 

Nyear 

Life  time  (years) 

10 

fmain 

Maintenance  factor  per  year 

0.1  (10%) 

lamor 

Interest  rate  per  year 

0.05  (5%) 

where  WE  is  the  work  recovered  by  the  expander.  It  is  through  the 
motor  power  that  the  transient  modeling  of  the  WRAS  is  accom¬ 
plished  via  both  a  transient  electrical  and  a  transient  torque  balance 
on  the  motor. 


Fig.  3.  Screw  compressor  performance  map  [8]. 


compressors,  they  are  a  better  choice  for  small  scale  fuel  cell 
systems.  The  compressor  work  rate  is  given  by 


Wc  = 


™C,airCp 


Tlcn 

Pc 


fp T 


((*-!)/*) 

-  1 


(4) 


2.2.3.  FPS 

Due  to  its  complex  structure,  the  FPS  model  requires  a  large 
amount  of  information  to  describe  its  kinetic  characteristics,  tran¬ 
sient  behavior,  and  design  procedures  which  cannot  be  repeated 
here  due  to  space  limitations.  Thus,  the  reader  is  referred  to  Kim 
[2]  for  more  details. 

3.  Economic  analysis 


where  rhc  air  is  the  mass  flow  rate,  cp  the  specific  heat,  T£n  and 
T°ut  the  inlet  and  outlet  temperatures,  p™  and  p£uf  the  inlet  and 
outlet  pressures,  rjc  the  compressor  efficiency,  and  k  the  specific 
heat  ratio.  Any  mechanical  efficiency  or  leakage  is  ignored  in  this 
model.  The  compressor  efficiency  varies  according  to  the  flow  rate 
and  the  pressure  ratio.  When  it  is  available,  dynamic  operation  of 
a  compressor  can  be  correctly  described. 

For  partial  loads,  the  map  in  Fig.  3  cannot  be  used  directly  for 
the  fuel  cell  system  because  this  map  is  for  a  250  kWe  system.  In 
order  to  describe  the  appropriate  operation  of  the  compressor  for 
the  5  kWe  PEMFC  system  developed  in  this  work,  the  map  is  scaled 
down  by  using  scaling  factors  for  the  mass  flow  rate  and  pressure 
ratio  as  follow: 

Nc=f(pcSFp,mcSF™)  (6) 

r]c=f(PcSFpc,  mcSF™ )  (7) 


Although  there  are  many  available  cost  models  for  a  variety  of 
energy  systems,  it  is  difficult  to  find  appropriate  ones  for  small 
scale  energy  systems  like  the  5kWe  PEMFC  system  considered 
here.  Moreover,  even  available  cost  models  for  small  scale  systems 
usually  show  significant  discrepancies  between  the  real  market 
prices  and  the  predicted  prices.  Using  unrealistic  cost  models  prej¬ 
udice  the  system  optimization  results.  Therefore,  it  is  necessary  to 
establish  and  use  appropriate  cost  functions  for  the  optimization 
process.  In  this  study,  cost  functions  for  the  5  kWe  PEMFC  system 
are  obtained  from  several  different  literature  sources  and  modified 
by  comparing  them  to  real  market  prices.  Table  1  summarizes  some 
assumptions  for  the  economic  analysis.  In  addition,  the  natural  gas 
price  is  obtained  from  the  Energy  Information  Administration  [4], 
and  an  average  residential  price  of  13.75$  left-3  is  used  for  evalu¬ 
ating  the  operating  costs.  Details  of  the  cost  functions  are  given  in 
Kim  [2]. 


where  SF™(=  pPD/p^esign)  is  the  scaling  factor  for  the  pressure  ratio 

and  SF™(=  mBcD /md^slgn)  that  for  the  mass  flow  rate.  pBD  and  mBD 
are  the  base  design  point  pressure  ratio  and  mass  flow  rate  factor, 
respectively.  pdesign  and  mdesign  are  the  design  pressure  ratio  fac¬ 
tor  and  mass  flow  rate  factor  determined  by  the  optimization.  The 
expander  and  DC  motor  designs  implemented  in  the  same  way  are 
described  in  detail  in  Kim  [2]. 

The  compressor  and  expander  are  connected  via  a  single  shaft 
so  that  the  extracted  power  of  the  expander  transfers  to  the  com¬ 
pressor.  The  DC  motor  and  compressor  are  linked  through  a  gear 
drive  as  shown  in  Fig.  4.  The  required  load  on  the  motor,  1%,  is 
given  by 

WM  =  WC-  WE  (8) 


Fig.  4.  Schematic  of  linkages  between  WRAS  components. 


4.  Optimization  strategy 

4.1.  Multi-level  optimization  strategy 

Conceptually,  the  decomposition  process  is  placed  in  between 
the  deterministic  model  and  the  optimizer  as  shown  in  Fig.  1. 
Decompositions  in  the  multi-level  approach  can  be  achieved  in 
four  ways,  i.e.,  disciplinary,  conceptual,  time,  and  physical  decom¬ 
position  [1].  In  this  study,  physical  decomposition  is  employed. 
Various  physical  decomposition  techniques  have  been  introduced 
in  the  literature  and,  in  general,  can  be  classified  as  methods 
of  either  local-global  optimization  (LGO)  or  iterative  local-global 
optimization  (ILGO).  A  dynamic  version  of  the  latter  also  exists  des¬ 
ignated  as  DILGO.  Both  ILGO  and  DILGO  have  been  developed  and 
their  effectiveness  validated  in  energy  system  synthesis/design  and 
operation/control  optimization  by  Munoz  and  von  Spakovsky  [5], 
and  Rancruel  and  von  Spakovsky  [6]. 

4.1.1.  DILGO  and  coupling  function  definition  for  the  dynamic 
PEMFC  system  optimization  under  uncertainty 

In  order  to  apply  the  DILGO  approach  to  the  dynamic  PEMFC 
system  optimization  problem,  the  system  is  decomposed  into 
sub-systems  connected  by  coupling  functions.  The  decomposition 
coincides  with  the  three  subsystems  described  above  (i.e.,  the  SS, 
FPS,  and  WRAS),  which  are  used  as  the  basis  for  replacing  the 
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Fig.  5.  Subsystem  boundaries  and  coupling  functions. 


system-level  optimization  problem  based  on  system  total  life  cycle 
cost  with  three  system-level,  unit-based  (SLUB)  optimization  sub¬ 
problems. 

The  system  configuration  in  Fig.  2  is  used  as  the  initial 
configuration  during  the  synthesis/design  and  operation/control 
optimization.  Each  subsystem  boundary  and  the  set  of  coupling 
functions  used  here,  which  represent  material/energy  flows  or 
required  operating  conditions  between  subsystems,  are  depicted 
in  Fig.  5.  Although  the  PES  is  not  modeled  here,  the  SS  and  the  PES 
are  coupled  by  a  controller  to  make  the  output  power  Pgr0ss  from  the 
SS  meet  the  power  requirement  Preq  of  the  load  demand.  Five  cou¬ 
pling  functions  are  defined  for  the  SLUB  optimizations:  the  anode 
and  cathode  side  pressures  p anode  and  P cathode  »  the  hydrogen  flow 
rate  hH2  from  the  FPS,  the  air  flow  rate  hAir  from  the  WRAS,  and  the 
motor  parasitic  power  EWRAS. 

The  SS  can  be  considered  as  the  main  subsystem  of  the  PEMFC 
system  because  the  FPS  and  the  WRAS  are  run  to  satisfy  the  SS 
requirements  (i.e.,  pressure,  hydrogen  and  air  flow  rates)  for  a  given 
load  demand.  To  generate  the  gross  electric  power  for  a  specific  load 
and  the  required  parasitic  power  EWRAS  >  the  SS  requires  a  certain 
amount  of  hydrogen  and  air  (i.e.,  hH2  and  hAir)  as  well  as  operat¬ 
ing  pressure  to  minimize  the  total  life  cycle  cost.  These  required 
hydrogen  and  air  flow  rates  and  pressure  become  the  base  pro¬ 
files  for  conducting  the  SLUB  optimizations  of  the  dynamic  FPS  and 
WRAS.  Even  though  stack  temperature  is  also  an  important  operat¬ 
ing  parameter,  inlet  temperatures  of  the  air  and  hydrogen  are  not 
considered  as  coupling  functions  because  these  temperatures  are 
managed  by  controllers  so  that  they  are  not  affected  by  the  other 
subsystems. 

To  evaluate  and  quantify  the  uncertainty  effects  on  system 
development,  three  uncertainty  factors  (i.e.,  with  respect  to  the 
load  profile,  the  cost  functions,  and  certain  inequality  constraints 
(e.g.,  that  on  the  CO  concentration))  are  considered  for  the  opti¬ 
mization.  For  example,  the  total  life  cycle  cost  can  be  expressed 
in  terms  of  expected  total  life  cycle  cost  with  a  certain  range,  i.e., 
it  can  be  expressed  by  deviation  of  the  expected  total  life  cycle 
cost  in  terms  of  the  standard  deviation.  The  RSA  method  has  been 
employed  to  quantify  the  uncertainties  in  the  simulations  since  it 
shows  high  fidelity  and  efficiency  for  the  large-scale  dynamic  sim¬ 
ulation  problems  [2,7]  without  the  large  computational  burdens 
imposed  by  traditional  stochastic  approaches  such  as,  for  example, 
Monte  Carlo  simulations. 

The  following  sections  provide  a  detailed  formulation  of  the 
SLUB  optimization  sub-problem  for  each  subsystem  with  consid¬ 
ered  uncertainties  in  the  objectives  and  constraints. 

4.1.2.  SLUB  optimization 

Total  life  cycle  cost  based  objective  functions  are  defined  for 
the  SLUB  optimization  sub-problems  and  summarized  in  Table  2. 
For  example,  the  vector  of  equality  constraints  7? FPS  represents  the 
thermodynamic,  kinetic,  and  geometric  models  of  the  FPS  and  the 
vector  of  inequality  constraints  C?  FPS  the  physical  limitations  on  the 


FPS  such  as  the  limit  on  CO  concentration.  Eq.  (9c)  indicates  that  the 
coupling  functions  hn2  and  panode  must  take  the  values  which  are 
found  by  solving  the  previous  SS  SLUB  optimization  sub-problem. 
Thus,  the  FPS  SLUB  optimization  sub-problem  is  solved  by  fixing 
the  inputs  and  optimizing  the  output,  i.e.,  in  this  case  the  natural 
gas  required. 

Eq.  (9)  is  the  objective  function  for  the  SLUB  optimization  of 
the  FPS.  The  SLUB  objective  function  C'FPS  consists  of  contributions 
from  the  FPS  (i.e.,  the  capital  cost  CFPS  (or  capital  investment  cost 
for  the  life  cycle)  of  the  FPS  and  the  cost  of  fuel  consumed  Coper), 
the  optimum  values  for  the  capital  costs  of  the  SS  C*s  and  WRAS 
^wras  obtained  from  solving  sub-problems  (11)  and  (12),  and  the 
variation  of  the  total  life  cycle  cost  in  terms  of  its  standard  devi¬ 
ation  (i.e.,  C™a[  UNC).  This  standard  deviation  is  expressed  by  Eqs. 
(9a)  and  (9b).  C^c  captures  uncertainties  in  the  SLUB  optimization 
of  the  FPS  based  on  uncertainties  in  the  cost  functions  (i.e.,  purchase 
cost)  of  the  FPS,  load  profile  (i.e.,  variations  in  the  hydrogen  require¬ 
ment),  and  fuel  cost.  C^*c,  C^J^S*  and  C^*c  represent  unit-level 
uncertainties  obtained  from  the  previous  SLUB  optimizations  of  the 
WRAS  and  SS  using  Eqs.  (lib)  and  (12b).  Operating  cost  defined  as 
Eq.  (10)  represents  total  hydrogen  requirement  which  is  described 
in  terms  of  cost.  In  order  to  estimate  the  life  cycle  operating  cost 
required,  the  hydrogen  consumption  rate  must  be  integrated  over 
the  entire  operating  time.  This  fuel  consumption  rate  is  described 
as  a  dynamic  load  profile  and  is  dealt  with  in  more  detail  in  Section 

4.2. 

As  to  the  inequality  constraint  for  the  CO  concentration  in  the 
FPS  SLUB  optimization,  GCo  is  the  mean  value  of  the  CO  concentra¬ 
tion  of  the  reformate  gas  after  the  PrOx  and  G ^  is  the  standard 
deviation  of  this  CO  concentration.  By  multiplying  the  standard 
deviation  by  2,  a  95%  confidence  interval  for  the  CO  concentration 
can  be  described.  In  addition,  the  parameters  u/oad,  Vfue[,  and  vCost 
are  the  variances  of  the  load  profile,  fuel  cost,  and  purchase  cost. 
An  8%  standard  deviation  of  the  load  profile  and  of  the  fuel  cost 
and  a  10%  standard  deviation  of  purchase  cost  are  assumed,  and  all 
probabilistic  distributions  of  the  uncertainty  sources  are  assumed 
as  normal  distributions. 

Finally,  a  sequential  scheme  for  the  DILGO  optimization 
approach  for  the  PEMFC  system  synthesis/design  and  opera¬ 
tion/control  under  uncertainty  is  proposed  as  shown  in  Fig.  6.  The 
DILGO  procedure  is  repeated  until  a  set  of  predefined  limits  (0.1% 
in  this  study)  on  the  SLUB  optimizations  are  reached.  Furthermore, 
note  that  the  optimal  control  architecture  and  its  control  gains  can 
be  determined  as  part  of  an  energy  system  synthesis/design  and 
operation/control  optimization  problem.  Wang  [8]  has  done  this 
for  the  optimal  control  architecture  design  of  this  PEMFC  system 
based  on  a  state-space  control  for  a  set  of  multi-input,  multi-output 
(MIMO)  controllers.  In  this  work,  the  controller  and  control  gains 
for  each  subsystem  are  obtained  from  Wang  [8]  where  the  gains  are 
obtained  from  a  set  of  preliminary  subsystem-level  optimizations 
by  Wang. 

4.2.  Load  profile 

In  order  to  conduct  the  optimization  of  the  proposed  dynamic 
PEMFC  system,  a  dynamic  load  profile  for  a  5  kWe  residential  build¬ 
ing  is  required.  For  this  purpose,  daily  electricity  consumption 
profiles  of  residential  buildings  in  2006  are  obtained  from  South¬ 
ern  California  Edison,  and  the  mean  and  standard  deviation  of  the 
load  profiles  are  determined.  In  this  study,  it  is  assumed  that  the 
power  demands  during  the  summer  and  winter  seasons  are  the 
most  critical  for  system  operation  because  more  power  is  required 
during  these  seasons  compared  to  the  spring  and  fall.  Therefore,  a 
two-day  load  profile  for  the  system  optimization  is  developed  in 
which  half  of  the  load  profile  represents  the  summer  and  the  other 


Table  2 

Objective  function  definition  of  the  SLUB  optimization  sub-problem  for  each  subsystem. 


Objective  function  (minimize) 


Constraints 


Uncertainty  terms 


FPS 

^FPS  =  ^FPS  +  Coper  +  C^aLUNC  +  Qs  +  ^WRAS 


(9) 


7?  FPS  —  0 

T^fps  <  0 


where  Coper 


nNGdt 


(10) 


"H2-n*2 
Pandoe  ~  P*anode 


=  0 


SS 


C'ss  -  ^ss  +  C^tal  UNC  +  C*ps  +  C^s  +  C*per  +  /  Xh2  AfiH2  dt 

Jt=  0 

C  C  C 

/  7pano(ie  Ap anode  dt  +  /  A,^jr  Afl^ir  dt  +  /  ^PCathode  ^P  cathode  dt 

J  t= 0  J  t= 0  J  t=  0 


(11) 


WRAS 

^  WRAS  =  CwRAS  +  C^UNC  +  C*s  +  C*ps  +  C*per 

fT  .  .  (12) 

+  /  A .t  AEu/ras  dt 

J  ewras  WKnil 

J  t= o 


Geo, 95%  =  Geo  +  2 GgJ  <10 


f?ss  =  0 
^ss  <  0 


[ Ewras  -  e'wras^  =  0 

WRAS  =  0 
WRAS  <  0 

"A ir-n*Air 

Pcathode  ^cathode 


r -FPS 

U  total,  UNC  ~~ 


V^C 


I  r^SS*  I  ("WRAS* 
UNC  ^  ^  UNC  ^  u  UNC 


(9a) 


CFPS 
!  INC 


rFPS+  _  rFPS- 
total, load  total, load 

2Sf^, 

load 


CFPS+  -  CFPS~  \ 

total, fuel  ^  total,  fuel  \  ^tpp$ 

2S™,  I  Vfi‘el 

fuel  J 


rFPS+  _  fFPS— 
total,  COST  total, COST 


2  8fps 

* UCOST 


2 


V 


FPS 

COST 


g£ 


( -FPS+  _  s-FPS— 

U  CO,  load  U  CO, load 

2  8FPS 

^ uload 


2 


V 


FPS 

load 


(9b) 


(9c) 


fSS 

U  total,  UNC  ~~ 


/  fSS  I  fFPS*  I  fWRAS* 
V  UNC  UNC  UNC 


(11a) 


rss  _ 
Lnwr  — 


rss+  _  rss- 

total,  COST  L total, COST 


(lib) 


fWRAS  _  /^WRAS~T7^s7~rf^ 
total,  UNC  ~  V  UNC  “r  UNC  “r  UNC 


(12a) 


'  (-WRAS+ 
fWRAS  _  I  ^  total, Air 
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25WRAS 
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1st  DILGO 

Iteration 


2nd  DILGO 

Iteration 


nth  DILGO 

Iteration 


Fig.  6.  Sequential  computational  scheme  for  the  DILGO  strategy. 


the  winter.  Fig.  7  describes  this  48-h  load  profile.  In  order  to  facil¬ 
itate  the  calculation  process,  the  peak  point  on  the  load  profile  is 
shifted  to  the  starting  point.  That  is  why  the  Oth  hour  and  24th 
hour  power  demands  are  different  in  Fig.  7.  As  seen  in  the  figure, 
the  mean  profile  represents  the  real  power  consumed  during  the 
summer  and  winter,  and  the  modified  profile  is  a  simplified  version 


Fig.  7.  Electricity  demand  profile  and  its  variance  over  48  h. 


of  the  mean  used  for  the  optimization  process.  Day-to-day  varia¬ 
tions  are  described  as  uncertainties  in  the  power  demand  of  the 
mean  load  profile.  Furthermore,  to  accommodate  operation  over 
a  whole  year,  the  modified  48-h  profile  is  multiplied  by  182.  Of 
course,  this  yearly  profile  does  not  result  in  one,  which  accurately 
reflects  the  off-seasons  of  spring  and  fall  since  it  simply  extends  the 
winter  and  summer  seasons  from  a  6-month  to  a  12-month  period. 
Nonetheless,  it  is  used  for  the  purpose  of  guiding  the  optimization 
towards  a  synthesis/design  more  conservative  than  might  occur 
if  the  spring  and  fall  seasons  had  been  employed  for  the  second 
six-month  period. 


5.  Results  and  discussion 

5.1.  Multi-level  optimization  results 

The  PEMFC  life  cycle  costs  and  the  uncertainties  on  these 
costs  for  each  iteration  of  the  DILGO  approach  are  presented  in 
Figs.  8  and  9.  Via  the  DILGO  process,  the  optimum  cost  of  each 
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Fig.  8.  Change  in  the  total  life  cycle  cost  of  the  PFMFC  system  during  the  DILGO 
procedure  with  a  95%  confidence  interval. 

subsystem  is  obtained  and  the  optimum  operating  cost  for  the 
entire  life  time  is  also  determined  via  the  dynamic  load  profile. 
Fig.  8  presents  the  total  life  cycle  cost,  while  Fig.  9  shows  the  opti¬ 
mum  costs  of  the  capital  for  each  subsystem  and  the  operating  cost 
for  each  iteration.  The  DILGO  procedure  continues  until  improve¬ 
ment  in  the  total  life  cycle  cost  is  below  0.1%.  The  global  optimum 
value  for  the  total  life  cycle  cost  of  the  PEMFC  system  obtained 
during  the  6th  iteration  of  the  DILGO  procedure  is  $38,077  with  a 
standard  deviation  of  $3932.  Thus,  the  optimum  total  life  cycle  cost 
of  the  PEMFC  system  is  expressed  as  $38,077  ±  $7864  with  a  95% 
confidence  interval.  Uncertainties  in  the  capital  cost  of  each  sub¬ 
system  and  the  operating  cost  are  presented  in  Fig.  9  using  error 
bars.  Preliminary  studies  [7]  indicated  that  the  mean  values  of  the 
optimum  total  life  cycle  cost  of  the  PEMFC  system  showed  little 
difference  whether  or  not  the  uncertainty  terms  were  included. 
However,  uncertainties  on  the  constraints  did  significantly  affect 
the  optimization  results.  In  particular,  uncertainty  considerations 
with  regard  to  the  CO  concentration  led  to  different  results  for  the 
SLUB  optimization  of  the  FPS.  That  effect  was  taken  into  account 
in  the  optimization  results  in  this  paper.  This  uncertainty  effect 


DILGO  iteration  number 

(a)  Capital  cost  of  the  SS 


O 

DILGO  iteration  number 

(c)  Capital  cost  of  the  WRAS 


on  optimization  is  explained  in  detail  by  the  authors  in  Kim  [2] 
and  Kim  et  al.  [7].  Readers  are  referred  to  those  references  for  a 
better  understanding  of  the  uncertainty  effects  on  PEMFC  system 
development. 

Significant  improvements  in  the  system-level  objective  function 
and  capital  costs  of  the  subsystems  are  observed  upon  completion 
of  the  second  DILGO  iteration,  and  there  are  no  significant  improve¬ 
ments  after  the  4th  iteration.  The  second  DILGO  iteration  predicts 
the  total  life  cycle  cost  of  the  PEMFC  system  to  be  $38,200  which  is  a 
1 8.44%  improvement  over  the  first  DILGO  iteration,  but  there  is  only 
a  0.26%  improvement  from  the  second  to  the  6th  DILGO  iteration. 
It  may  appear  that  the  optimization  of  the  PEMFC  is  already  com¬ 
pleted  after  the  second  or  the  third  DILGO  iteration  because  there 
is  little  improvement  after  either.  However,  to  determine  whether 
or  not  the  system-level  optimization  using  DILGO  has  completed, 
the  stabilization  of  all  the  coupling  functions  must  be  evaluated  as 
well. 

Fig.  10  depicts  this  stabilization.  Fig.  10a  and  b  presents  the 
required  hydrogen  and  air  flow  rate  at  the  full  load  condition.  Each 
point  in  the  graphs  coincides  with  the  starting  point  of  the  dynamic 
profile  in  Fig.  7.  The  required  hydrogen  flow  rate  and  the  required 
air  flow  rate  stabilize  at  0.044  mol  s-1  and  0.208  mol  s-1  at  the  full 
load  condition,  respectively.  Another  coupling  function,  the  motor 
parasitic  power  seen  in  Fig.  10c,  oscillates  with  each  DILGO  itera¬ 
tion  and  begins  to  stabilize  at  around  420  W  after  the  6th  DILGO 
iteration.  As  to  the  stack  operating  pressure,  it  is  one  of  the  cou¬ 
pling  functions  as  well  as  one  of  the  operation  decision  variables. 
The  pressure  plays  an  important  role  in  PEMFC  system  operation 
because  it  affects  PEMFC  system  performance  and  also  determines 
the  operating  pressure  of  the  FPS  and  WRAS.  A  final  optimum  oper¬ 
ating  stack  pressure  appears  to  stabilize  at  2.22  bar  in  the  6th  DILGO 
iteration  as  seen  in  Fig.  1  Od.  Thus,  it  can  be  inferred  from  the  figures 
that  stabilized  coupling  functions  may  be  obtained  after  one  or  two 
more  DILGO  iterations. 

The  global  optimum  synthesis/design  and  operation  decision 
variable  values  of  the  PEMFC  system  appear  in  Table  3.  The  FPS  SLUB 
optimization  eliminates  the  first  pre-heater,  P_HX1,  and  the  third 
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Fig.  10.  Changing  coupling  functions  during  the  DILGO  procedure. 


heat  exchanger,  HX3  from  the  initial  configuration  in  Fig.  2.  The 
optimal  design  of  the  SS  predicts  0.127  m2  for  the  single  cell  active 
area  (3.7  m2  of  total  cell  active  area),  while  the  initial  design  of  the 
SS  has  a  single  active  area  of  0.048  m2  (2.0  m2  total  cell  active  area). 
The  PEMFC  stack  efficiency  increases  with  increasing  stack  (active 
area)  size  because  the  stack  can  operate  at  lower  current  densities 
as  seen  in  Fig.  1 1.  In  this  figure,  the  cell  current  density  and  hydro¬ 
gen  consumption  rate  decrease  with  increasing  cell  active  area.  Less 
cell  current  density  produces  higher  cell  voltage  to  generate  the 
same  amount  of  power.  Thus,  the  larger  size  fuel  cell  stack  con¬ 
sumes  less  hydrogen.  This  in  turn  reduces  the  operating  cost  of  the 
PEMFC  system.  However,  as  the  fuel  cell  stack  size  increases,  the 
capital  cost  of  the  SS  increases  as  well.  Therefore,  a  compromise 
between  the  operating  cost  of  the  PEMFC  system  and  the  capital 
cost  of  the  SS  with  regard  to  the  stack  size  must  be  determined 
by  the  system-level  optimization  process.  As  seen  in  Table  3,  this 
compromise  results  in  a  cell  active  area  of  0.127  m2  with  the  opti¬ 
mum  number  of  cells  being  29.  Even  though  the  larger  cell  size  may 
cause  some  technical  difficulties  in  terms  of  stack  assembly,  ther¬ 
mal  management,  and  fuel  distribution  in  the  cell,  these  are  not 
considered  in  this  work  because  there  is  no  available  information 


Fig.  11.  Relationships  between  the  single  cell  active  area  and  current  density  and 
the  H2  consumption  rate  at  full  load. 


about  size  limitations  on  stack  design  in  the  literature.  If  it  becomes 
available,  it  must  be  taken  into  account  as  a  constraint. 

The  optimum  purchase  cost  of  each  subsystem  and  its  cost 
breakdown  are  presented  in  Table  4  based  on  a  production  rate  of 
500,000  manufactured  units  per  year.  Note  that  these  costs  are  pro¬ 
jected  future  costs  which  reflect  a  mature  technology  in  which  the 
major  technical  and  lifetime  issues  have  been  largely  resolved  and 
the  volume  of  units  manufactured  has  reached  levels  sufficiently 
large  to  drop  the  cost  of  this  technology  to  levels  which  make  it 
competitive  with  more  traditional  technologies.  With  this  qualifi¬ 
cation,  the  optimum  purchase  cost  of  the  PEMFC  system  is  $3720, 
and  the  optimum  purchase  costs  of  the  FPS,  SS,  and  WRAS  make 
up  51%,  39%,  and  10%  of  the  system  purchase  cost,  respectively.  As 
shown  in  the  table,  the  cost  of  the  reactors  (the  SMR,  HTSR,  LSTR, 
and  PrOx  reactor)  is  around  62%  of  the  optimum  purchase  cost  of  the 
FPS,  that  of  the  three  heat  exchangers  (P_HX2,  HX1,  and  HX2)  19%, 
and  that  of  the  auxiliary  parts  (e.g.,  valves,  pipes,  controllers,  etc.) 
1 7%.  The  remaining  2%  is  the  steam  generator.  The  stack  purchase 
cost  takes  82%  of  the  optimum  purchase  cost  of  the  SS,  and  89% 
of  that  is  for  purchasing  the  electrodes  which  include  the  anode, 
cathode,  and  gas-diffusion-layers.  The  purchase  cost  of  the  auxil¬ 
iary  parts  such  as  pipes,  valves,  and  controllers  for  the  SS  represents 
8%  of  the  SS  purchase  cost.  It  should  be  noted  that  some  of  auxil¬ 
iary  parts  for  the  SS  are  also  common  with  the  WRAS  because  both 
subsystems  are  coupled  with  each  other  via  these  valves  and  con¬ 
trollers.  Thus,  the  purchase  cost  of  these  parts  for  the  SS  and  WRAS 
is  not  considered  in  the  purchase  cost  of  the  WRAS  but  in  that  of 
the  SS. 


5.2.  Dynamic  response  of  the  SS 

The  dynamic  responses  of  the  cell  current  density  and  cell  volt¬ 
age  over  the  48  h  of  operation  are  depicted  in  Fig.  12.  In  this  figure, 
the  optimum  responses  are  compared  to  the  non-optimum  (i.e.,  ini¬ 
tial  synthesis/design  and  operating  conditions)  responses.  As  seen 
in  the  figure,  the  optimum  cell  current  density  is  lower  than  the 
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Table  3 

Optimum  synthesis/design  and  operation  decision  variable  values. 


Synthesis/design  decision  variables 

Component 

Decision  variables 

Initial 

Optimum 

Component 

Decision  variables 

Initial 

Optimum 

value 

value 

value 

value 

FPS 

SMR 

LSMR  (reactor  length,  m) 

1.0 

1.15 

P_HX2 

Lp*HX2  (heat  changer  plate  length,  m) 

0.4 

0.300 

d™R  (reactor  tube  diameter,  m) 

0.02 

0.0237 

2  (channel  height,  m) 

0.004 

0.0038 

HTSR 

IHtsr  (reactor  length,  m) 

0.5 

0.590 

Npjjx2  (Number  of  plates  for  each  side) 

4 

4 

dj1™  (reactor  tube  diameter,  m) 

0.09 

0.0663 

HX1 

(heat  changer  plate  length,  m) 

0.25 

0.540 

LTSR 

ILtsr  (reactor  length,  m) 

0.6 

0.680 

(channel  height,  m) 

0.004 

0.0039 

dLTSR  (reactor  tube  diameter,  m) 

0.09 

0.0708 

(number  of  plates  for  each  side) 

3 

4 

PrOx 

LPr0x  (reactor  length,  m) 

0.5 

0.540 

HX2 

l{J*2  (heat  changer  plate  length,  m) 

0.3 

0.197 

d^0x  (reactor  tube  diameter,  m) 

0.08 

0.0670 

M!x2  (channel  height,  m) 

0.004 

0.0039 

SG 

LSG  (steam  generator  tube  length,  m) 

1.0 

1.71 

N™2  (number  of  plates  for  each  side) 

3 

4 

P_HX1* 

Lp^x i  (heat  changer  plate  length,  m) 

0.2 

N/A 

HX3* 

3  (heat  changer  plate  length,  m) 

0.3 

N/A 

hp^Hx i  (channel  height,  m) 

0.004 

N/A 

3  (channel  height,  m) 

0.004 

N/A 

Np^nx i  (number  of  plates  for  each  side) 

4 

N/A 

N™3  (number  of  plates  for  each  side) 

2 

N/A 

WRAS 

Motor 

^des'Z  (desi§n  torque,  Nm) 

3 

2.10 

Expander 

PdeZnder  (design  pressure  ratio) 

2.0 

2.17 

rPmd£gn  (design  rotational  speed,  rpm) 

2400 

1457 

mfp™der  (design  mass  flow  factor) 

0.1 

0.076 

Compressor 

PCd°e™gnSSOr  (desi§n  Pressure  ratio) 

2.2 

2.18 

Unit: 

m Compressor  &  ^ Expander  (|<gT0.5 /s  bar)  JnMotor 
design  ",  design  1  a)'  ^design 

SS 

Stack 

m^T  (design  mass  flow  factor) 

0.1 

0.083 

(Nm) 

Lss  (cell  length,  m) 

0.218 

0.357 

Stack 

Nss  (number  of  cells) 

42 

29 

Operation  decision  variables 

Variable 

Initial  value 

Optimum  value 

Variable  Initial  value 

Optimum  value 

£s/c  (steam- to-carbon  ratio) 

3.5 

2.92 

tHX  3 
$AIR 

(cooling  air  flow  ratio  of  HX3)  0.12 

N/A 

£ch4  (fuel  feed  ratio  between  SMR  and  combustor) 

0.42 

0.305 

hPrOx 

^AIR 

(cooling  air  flow  ratio  of  PrOx)  0.40 

0.300 

%airR  (cooling  aif  flow  ratio  of  HTSR) 

0.13 

0.186 

T™R  (inlet  gas  temperature  of  HTSR,  I<)  610 

612 

%air2  (cooling  air  flow  ratio  °f  HX2) 

0.28 

0.365 

pss  ( stack  operating  pressure,  bar)  1 .8 

2.2 

(cooling  air  flow  ratio  of  LTSR) 

0.07 

0.149 

*  Component  has  been  eliminated  by  the  optimization 

Table  4 

Optimum  purchase  cost  based  on  500,000  units/year. 


SS($) 

Stack 

Membrane 

81.1 

Bipolar  Plate 

46.2 

Electrode 

1046.7 

Assembly 

31.5 

Stack  purchase 

1205.5 

Humidifier 

57.5 

Cooling  cycle 

76.6 

Auxiliary  parts 

114.8 

SS  purchase 

1454.4 

FPS ($) 

Reactors 

1164.9 

Steam  Generator 

41.4 

Heat  exchangers 

365.7 

Auxiliary  parts 

314.4 

FPS  purchase 

1886.5 

WRAS  ($) 

Compressor 

225.6 

Expander 

108.5 

Motor 

45.1 

WRAS  purchase 

379.2 

Time  (s) 


Fig.  12.  Dynamic  responses  of  the  cell  current  density  and  cell  voltage  of  the  opti¬ 
mum  and  initial  syntheses/designs. 


non-optimum  cell  current  density.  Thus,  the  optimum  cell  voltage 
is  higher  than  the  non-optimum  cell  voltage  to  generate  the  same 
amount  of  electrical  power.  As  explained  in  the  previous  section, 
the  larger  size  fuel  cell  consumes  less  hydrogen,  so  operating  costs 
are  reduced.  Through  the  system-level  optimization,  about  47%  of 
the  initial  cell  current  density  is  predicted  as  the  optimum  cell  cur¬ 
rent  density.  The  optimum  dynamic  cell  current  density  at  the  full 
load  condition  (at  Os)  is  0.1 8  A  cm-2  as  compared  to  0.38  A  cm-2 
for  the  initial  synthesis/design.  At  40%  of  full  load  (at  43,200  s),  it 
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Fig.  13.  Dynamic  responses  of  the  compressor,  expander,  and  motor  power  for  the 
optimum  and  initial  designs. 
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Fig.  15.  Comparison  of  the  optimum  and  initial  (non-optimized)  FPS  efficiencies. 

5.4.  Optimum  system  efficiency 


- Efficiency  (optimum) 

-  -  -  -  Efficiency  (initial) 


is  about  0.07  A  cm-2  as  opposed  to  the  about  0.1 4  A  cm-2  for  the 
initial  system. 

5.3.  Optimum  design  and  dynamic  responses  of  the  WRAS 

The  optimum  dynamic  responses  of  the  compressor,  expander, 
and  motor  work  rate  are  compared  in  Fig.  13  with  the  dynamic 
responses  of  the  initial  design  over  the  48  h  of  operation.  The  initial 
non-optimized  design  of  the  compressor  requires  95  W  more  of 
power  at  the  full  load  condition  (at  Os)  and  40 W  more  at  40%  of 
full  load  (at  43,200)  than  the  optimum  compressor.  In  the  case  of 
the  expander  unit,  the  initial  non-optimized  design  of  the  expander 
recovers  35  W  more  at  full  load  and  9W  more  at  40%  of  full  load 
than  the  optimum  design.  Therefore,  taking  into  account  the  initial 
non-optimized  motor  design,  the  initial  non-optimized  design  of 
the  compressor  needs  an  additional  61 W  of  work  at  full  load  and 
42  W  more  at  40%  of  full  load.  This  is  the  case  because  for  the  initial 
motor,  which  is  non-optimized,  the  additional  motor  work  rate  is 
slightly  higher  than  the  difference  between  the  additional  work  for 
the  compressor  (i.e.,  95  W  at  full  load  and  40  W  at  40%  of  full  load) 
and  the  additional  work  recovered  by  expander  (i.e.,  35  W  at  full 
load  and  9  W  at  40%  of  full  load). 

The  performance  of  the  WRAS  can  also  be  expressed  in  terms 
of  the  work  recovery  rate  which  is  the  rate  of  the  work  recovered 
by  the  expander  for  the  required  work  of  the  compressor.  This  is 
described  in  Fig.  14  which  compares  the  work  recovery  rate  of  the 
optimum  design  and  that  of  the  initial  non-optimized  design.  The 
optimum  WRAS  recovers  about  41%  with  the  required  compressor 
work  by  the  expander  unit  at  full  load,  and  about  28.5%  at  40%  of 
full  load.  The  work  recovery  rate  of  the  optimum  WRAS  is  from  0.5 
point  percent  (at  full  load)  to  1.5  point  percent  (at  40%  of  full  load) 
higher  than  that  of  the  initial  non-optimized  WRAS  as  seen  in  the 
figure. 


Fig.  14.  Comparison  of  the  optimum  and  initial  work  recovery  ratio  by  the  expander. 


The  efficiencies  of  the  FPS,  r]FPS,  and  the  PEMFC,  Ppemfc,  system 
are  evaluated  in  this  section  and  are  defined  as  follow: 


rhH2LHVH2 

(m^  +  <H4)LWcH4+<2mVH2 


(13) 


P  PEMFC  = 


Enet 


(14) 


where  rh^,  and  are  the  methane  flow  rates  to  the  SMR 
and  burner  and  the  hydrogen  flow  rate  to  the  burner,  respectively. 
Un-reacted  hydrogen  in  the  stack  is  fed  into  the  burner.  As  a  result 
of  this,  the  overall  system  efficiency  is  improved.  LHVC h4  and  LHV h2 
are  the  lower  heating  values  of  the  methane  and  hydrogen,  and  Enet 
is  the  system  net  power. 

Fig.  1 5  shows  the  dynamic  optimum  FPS  efficiency  over  the  48  h 
of  operation  and  compares  it  with  that  for  the  initial  non-optimized 
FPS.  The  optimum  FPS  efficiency  ranges  from  about  85.6  to  86.1% 
through  the  entire  operating  regime,  while  that  for  the  initial  non- 
optimized  FPS  varies  from  about  77.9  to  78.5%.  Initially,  a  0.42  fuel 
ratio  is  assumed  and  0.305  is  obtained  from  the  optimization,  which 
principally  leads  to  an  8%  efficiency  enhancement  of  the  FPS. 

The  optimum  overall  system  efficiency  of  the  PEMFC  system 
shown  in  Fig.  16  is  compared  with  that  of  the  initial  non-optimized 
PEMFC  system.  The  optimum  efficiency  of  the  PEMFC  system 
remains  fairly  steady  at  46.5%  throughout  the  operating  regime, 
while  that  for  the  initial  non-optimized  system  ranges  from  about 
36  to  36.5%.  Thus,  the  overall  system  efficiency  increases  by  1 0  per¬ 
centage  points  through  the  optimization  process.  The  main  reason 
for  the  improvement  is  that  the  hydrogen  fuel  consumption  of  the 
optimum  PEMFC  system  is  reduced  significantly  compared  to  that 
for  the  initial  non-optimized  PEMFC  system.  The  initial  design  of 
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Fig.  16.  Comparison  of  the  optimum  and  initial  (non-optimized)  PEMFC  system 
efficiencies  with/without  the  expander  unit. 
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the  stack  requires  0.0496  mol  s-1  of  hydrogen  at  full  load  but  only 
0.0436  mol  s-1  for  the  optimum  stack. 

The  effects  of  the  work  recovery  unit  using  the  expander  are 
evaluated  in  terms  of  the  PEMFC  system  efficiency  in  Fig.  16  as 
well.  As  seen  in  the  figure,  the  system  efficiency  is  enhanced  by  3.4 
percentage  points  by  using  the  work  recovery  unit  for  the  optimum 
PEMFC  system,  and  by  only  1.2  percentage  points  for  the  initial 
non-optimized  PEMFC  system. 

6.  Conclusions 

Decomposition  techniques  are  a  very  useful  strategy  for  large- 
scale  energy  system  optimization.  In  particular,  the  overall  system 
model  with  its  complex  control  architecture  and  very  detailed 
component  models  results  in  a  large  computational  burden  and 
quite  often  simulation  failure  during  a  dynamic  optimization. 
Fortunately,  the  optimization  of  the  complex  system  model  and 
control  architecture  can  be  treated  as  a  set  of  approximate 
smaller  unit  (subsystem)-level  optimization  problems  via  a  decom¬ 
position  strategy  such  as  DILGO.  This  physical  decomposition 
approach  is  successfully  applied  to  the  synthesis/design  and  oper¬ 
ation/control  optimization  of  the  dynamic  PEMFC  system,  and 
the  global  optimum  solution  is  found  within  6  iterations  by 
DILGO. 

The  optimum  synthesis/design  of  the  PEMFC  system  developed 
here  shows  a  fairly  steady  and  remarkably  high  overall  system  effi¬ 
ciency  at  46.5%  throughout  the  operating  regime  (i.e.,  full  load  to 
40%  of  full  load).  The  reason  is  that  all  design  issues  (i.e.,  synthe¬ 
sis/design  and  operation/control)  are  taken  into  account  in  one 
optimization  problem  so  that  all  the  optimum  values  are  achiev¬ 
able.  Moreover,  system  efficiency  can  be  significantly  enhanced  by 
using  the  work  recovery  unit  for  the  PEMFC  system,  and  this  study 
shows  a  3.4  percentage  points  increase  in  system  efficiency. 

As  to  the  operating  cost,  it  dominates  the  life  cycle  cost  of  the 
PEMFC  system  during  the  optimizations.  Thus,  minimizing  the  fuel 
consumption  rate  of  the  system  is  of  great  importance.  The  opti¬ 
mum  synthesis/design  of  the  FPS  plays  a  very  important  role  in  the 
entire  PEMFC  system  synthesis/design  and  operation/control  opti¬ 
mization  problem  because  most  of  the  additional  fuel  consumption 
in  the  PEMFC  system  occurs  during  the  fuel  processing  (e.g.,  the 
hydrogen  oxidation  reactions  in  the  preferential  oxidation  (PrOx) 


reactor,  fuel  for  the  combustor,  poor  reforming  performance,  etc.). 
The  large  cell  active  area  also  favors  the  reduction  in  the  fuel  con¬ 
sumption  rate. 

Finally,  the  approach  presented  here  is  an  effective  way  of 
dealing  in  a  single  problem  with  dynamic  fuel  cell  system  synthe¬ 
sis/design  and  operation/control  optimization  under  uncertainty. 
It  is  a  general  approach  which  can  be  applied  to  the  optimal  devel¬ 
opment  of  any  energy  system. 
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